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ABSTRACT 

Combing the three-dimensional radiative transfer (RT) calculation and cosmo- 
logical SPH simulations, we study the escape fraction of ionizing photons (/esc) of 
high-redshift galaxies at z ~ 3 — 6. Our simulations cover the halo mass range of 
Mh = 10^ - 10^2^^^ postprocess several hundred simulated galaxies with the 
Authentic Radiative Transfer (ART) code to study the halo mass dependence of /osc- 
In this paper, we restrict ourselves to the transfer of stellar radiation from local stellar 
population in each dark matter halo. We find that the average /esc steeply decreases 
as the halo mass increases, with a large scatter for the lower mass haloes. The low 
mass haloes with Mh ^ IO^Mq have large values of /osc (with an average of ^ 0.4), 
whereas the massive haloes with Mh ^ lO^^Af0 show small values of /osc (with an 
average of ^ 0.07). This is because in our simulations, the massive haloes show more 
clumpy structure in gas distribution, and star-forming regions are embedded inside 
these clumps, making it more difficult for the ionizing photons to escape. On the 
other hand, in low mass haloes, there are often conical regions of highly ionized gas 
due to the shifted location of young star clusters from the center of dark matter halo, 
which allows the ionizing photons to escape more easily than in the high-mass haloes. 
By counting the number of escaped ionizing photons, we show that the star-forming 
galaxies can ionize the intergalactic medium at z = 3 — 6. The main contributor to 
the ionizing photons is the haloes with Mh < 10^° owing to their high /osc- The 
large dispersion in /osc suggests that there may be various sizes of Hii bubbles around 
the haloes even with the same mass in the early stages of reionization. We also ex- 
amine the effect of UV background radiation field on /osc using simple, four different 
treatment of UV background. 

Key words: radiative transfer - ISM: dust, extinction - galaxies: evolution - galaxies: 
formation ~ galaxies: high-redshift - methods: numerical 



1 INTRODUCTION 

Observations of cosmic microwave background radiation 
provides a we alth of informatio n on the cosmic reioni zation 
history ( e.g.. |Page et al" ' '2007'; iDunklev et al] 120091 '). For 
example, IKomatsu et al. (201C1) showed that the reion- 
ization occurred at z ~ 10.5 assuming an instantaneous 
reionization scenario. However, the detailed history of 
reionization and the nature of ionizing sources are not 
yet fully understood. Since the UV background (UVB) 
radiation can heat up the interstellar medium (ISM) 
to ~ 10^ K and disturb star formation, UVB coupled 
with the ionization history of the universe significantly 
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influe n ces the galaxy f ormati o n (e.g., Susa fc Umemural 
2000l: lUmemura et all I2001I: ISusa fc Umemural 12004 



Okamoto et al.ll2008l : iHasegawa et al.ll2009l ). Therefore it is 



very important to study the UVB intensity and the nature 
of ionizing sources. 



iHaardt fc Madaul ([1991) pointed out that the UVB is 
dom inated by quasar s at z < 4. Using the SDSS sam- 



ple, iFan et all ^00^ showed that the bright-end slope 
of the quasar luminosity function at 2 >, 4 are consider- 
ably steeper than that at lower redshifts, and concluded 
that the quasars cannot maintain the ionization of IGM 
at z >j 4. Subsequently, much argument have been fo- 
cused on the possibility that the IGM is ionized mainly 
by the UV radiation from high-redshift (hereafter high-z) 
star-forming galaxies (e.g.. iFan et al.|[2006l : iBouwens et al.l 
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I2OO7I : iGnedinllioog ). The key quantity in determining the 
IGM ioni zation rate is the escape fraction of ionizing pho- 
tons (e.g.. iRazoumov fc Sommer^ Larsenll2006l :[ Gnedin et al.l 
I2OO8I ). which is the number ratio of photons escaping from 
a galaxy to the intrinsically radiated photons by stars. This 
parameter controls the contribution to the UVB intensity 
from star-forming galaxies. In this work, we examine the 
values of /esc in high- 2; star-forming galaxies. 

The re are severa l obser vational constraints on f^Bc at 
2 ~ 3. ISteidel et all (|200ll ) found /osc.roi > 0.5 from the 
composite spectrum of 29 Lyman Break Galaxies (LBGs) 
at 2 ~ 3, where /eac.rei is the relative fraction of escaping 
Lyman continuum (900 A) photons relative to the fraction 
of escaping non-ionizing UV (1500 A) photons. It is usually 
defined as 



esc,rel 



(L1500/L900)int 
(F1500/F900)obs 



exp(r9oo ) 



(1) 



where (F1500/F900)obs, (L1500/L900)int and r^g^ repre- 
sent the observed 1500 A/900 A flux density ratio, the in- 
trinsic 1500 A/900 A luminosity density ratio, and the line- 
of-sight opacity of the IGM for 900 A photons, respectively. 
Equation ([T]) compares the observed flux density ratio (cor- 
rected for the IGM opacity) with the models of UV spectral 
ene rgy distribution o f star- formin g galaxies. 

iGiallongo et all \2mi ) and llnoue et al.l (|2005l ') esti- 
mated t he upper limit o f /eac .rei ^ 0.1 — 0.4 for some LBGs 
at 2 3. IShaplev et al.l l)2006l ) directly detected the escaped 
ionizing photons from 2 LBGs in the SSA22 field at z = 3.1, 
and estimat e d the average value of /esc.rei = 0.14. Moreover, 
llwata et al.] (|2009l l successfully detected the Lyman contin- 
uum emission from 10 Ly-a emitters (LAEs) and 7 LBGs 
within a sample of 198 LAEs and LBGs in the SSA22 field. 
They showed that the mean value of /esc.rei for the 7 LBGs 
is 0.11 after correcting for dust extinction, and 0.20 if the 
IGM absorption is taken into account. 

In the early theoretical works, some authors stud- 
ied the feac with id eally modelled galaxies. For example, 
iDove fc Shuiil (|l994l ') estimated the /esc of Milky Way type 
galax y using a semi-analytic method, and reported /oac ~ 
0.07. iRicotti fc ShulJ (|200(]| ) investigated the dependence 
of /esc on various physical quantities, such as the collapse 
redshift and star formation e fficie ncy using a semi-an alytic 
method. IWood fc Loebl (|2000l ) and lCiardi etakl (|2002l ) stud- 
ied the effect of inhomogeneous structure of gas on /eac, and 
showed that /esc increases in clumpy systems by a factor 
of > 2 than in a homogeneous gas distribution. iDove et al.l 
(|200(j ) investigated the influence of bubbles made by super- 
novae on /es c using a sem i -analy tic method. Using numerical 
simulations, iFuiita et al.1 (|2003l ) studied the effect of super- 
novae feedback, and reported a high /esc (> 0.2) for a disk 
galaxy with Mn = 10** - IO^Mq. 

Theoretical studies in a more fully cosmological envi- 
ronment can be performed by combining cosmological hy- 
drodynamic simulations of galaxy formation and a three- 
dimensional r adiat ive transfer calculation. For example, 
lYaiima et al.l (|2009l . hereafter Y0 9) post-processed the eule - 
rian hydrodynamic simulation of lMori fc Umemural (|2006l ) 
with RT, and showed that the galaxies in an isolated halo 
of Mh = 10"Mo can have relatively large values of /eac = 
0.17 — 0.47. Moreover they found that /esc decreases grad- 



ually as a function of time owing to the dust pollution and 

the shifting star formation sites^ 

On the other hand, IRazoumov fc Sommer-LarsenI 

120061 ) examined the escape fractions of two galaxies in a 
cosmological SPH simulation from z = 3.8 to 2.4, which 
later become Milky Way type disk galaxies at z = 0. They 
found small values of /eac < 0.1, in disagreement with Y09. 
However they also reported that /esc decreases with redshift 
fro m z = 3.8 to 2.4, in qualitative agreem ent with Y09. 

IRazoumov fc Sommer-La^ (2OIOI ) further examined 
the /esc of star-forming galaxies in a wide mass range (M^ = 
lO'^'^-lO^ '^M©) at z = 4-10, and found that /esc decreases 
steeply as the halo mass increases in th eir cosmologic a l SPH 
sim ulations, in contrast to the work bv lGnedTn" et al.1 (|2008h 
and lWise fc Gen! (|2009l l. 

Using cosmological AMR simulations, I Gnedin et al.l 

(|2008l ) reported that haloes with Mh = lO" - IO^^Mq 
have /eac = 0.01 — 0.03, and much lower /esc for lower mass 
haloes with Mh = 10^" - 10^^ Mq. Their results suggest 
that /esc increases with halo mass, at least in the range of 
Mh = 10'° - 10"Mfl. 

IWise fc Gen! (|2009l ) extracted 10 haloes with masses 
Mh = 3 X 10*^ - 3 X W^Mq at 2 = 8 from cosmological 
AMR radiation hydrodynamic simulations, and examined 
the escape fraction of ionizing photons. They found that 
/esc fluctuates rapidly on a time-scale of a few to 10 Myrs 
depending on the star formation rates, and varies widely 
from almost zero to nearly unity. They found /esc ~ 0.4 for 
a normal IMF for the haloes with Mh = W'^'^ - 10^'^Mq, 
but /esc ~ 0.05 — 0.1 for lower mass haloes, disregarding the 
effect of dust. 

Although the halo mass dependence of /oac is very im- 
portant for the study of cosmic reionization, there are signif- 
icant differences in the theoretical estimates from cosmologi- 
cal hydrodynamic simulations as described above. In partic- 
ular, in these previous works, the number of studied haloes 
has been very small 10), therefore it has been difficult 
to gauge the halo mass dependence of /eac- In the present 
paper, we calculate the values of /eac for a much larger num- 
ber of haloes (several hundreds) in cosmological volumes of 
comoving 10— 100 /i^^Mpc, and examine its halo mass de- 
pendence. In addition, we study the effects of interstellar 
dust and UVB radiation on /eac. 

The outline of the paper is as follows. In §(21 the models 
and numerical methods are described. We present the results 
on escape fractions in § [31 and discuss the dust effect and 
the contribution of star-forming galaxies to the reionization 
of the universe in § |4] We then summarise in § (5] 



2 MODEL AND METHOD 
2.1 Simulations 

We use an updated and modified version of the Tree- 
particle-mesh (TreePM) smoothed particle hydro dynamics 
(SPH ) code GADGET-3 (originally described in ISpringell 
l2005h . The SPH calculation is performed based on the en - 
tropy conservative formulation (jSpringel fc Hernauistl2()o3 ') . 
Our fid ucial code includes radiative cooling by H, He, and 
metals (|Choi fc Nagaminl l2009bh . star formation, super- 
nova feedback, a phenomenological model for galactic winds 



© 2008 RAS, MNRAS 000.[HfT3l 



Escape of ionizing photons 3 



Series 


Box-size 






m-DM 


m 


gas 


e 




N144L10 


10.00 


2 


X 144^ 


1.97 X 10^ 


4.04 


X 


108 


2.78 


2.75 


N216L10 


10.0 


2 


X 216^ 


5.96 X 10^ 


1.21 


X 


10'^ 


1.85 


2.75 


N400L100 


100.0 


2 


X 400^ 


9.12 X 10® 


1.91 


X 


10® 


6.45 


0.0 



Table 1. Series of simulations employed for the present study. The box-size is given in units of /i~^Mpc, Np is the particle number of 
dark matter and gas (hence X 2), ttidm and rrtgas are the masses of dark matter and gas particles in units of h~^MQ, respectively, e 
is the comoving gravitational softening length in units of h^^kpc, and is the ending rcdshift of the simulation. The value of e is a 
measure of spatial resolution. 



and a sub-resolution model of 
multiphase ISM (jSpringel fc HernauistllioO^ ). We also in- 
clude the heating by a uniform UVB, which we will discuss 
more in ^ 12.41 



In this multiphase ISM model, high-density ISM is pic- 
tured to be a two-phase fluid consisting of cold clouds in 
pressure equilibrium with a hot ambient phase. Cold clouds 
grow by radiative cooling out of the hot medium, and this 
material forms the reservoir of baryons available for star for- 
mation. The star formation rate (SFR) is estimated for each 
gas particle that have densities above the threshold density, 
and the star particles are spawned statistically based on the 
SFR. For th e star formation model, th e "Pressure model" 
described in IChoi fc Nagamind (|2009al ) is being used. This 
model estimates the SFR based on the local gas pressure 
rather than the gas density, and implicitly considers the ef- 
fect of H2 formation. 

The simulations used in this paper also uses the new 
" Multicomponent Variable Velocity" wind model developed 
bv lChoi fc Nagamind (|2010l '). This new wind model is based 
on both the energy-driv en wind and the mom entum-driven 
wind model discussed bv lMurrav et al.| (|2005l ). and the wind 
speed in this model depends on the galaxy stellar mass 
and SFR. It gives more favorable results, and agrees bet- 
ter with the observations of, e.g., cosmic Civ mass den- 
sity and IGM temperature than the previous model with 
a constant wind s peed. To enable this new wind model, 
IChoi fc Nagamind (|201(]| ') implemented a on-the-fiy group- 
finder into GADGET-3 to compute the galaxy masses and 
SFRs while the simulation is running. The group-finder, 
which is a si mplified varian t of th e SUBFIND algorithm 
developed by ISpringel et al] (|200ll '). identifies the isolated 
groups of star and gas particles (i.e., galaxies) based on the 
baryonic density field. T he detailed procedure o f this galaxy 
grouping is described in iNagamine et all (|2004l ). The outer 
baryonic density threshold for a galaxy is O.Olpth, where pth 
is the threshold density for star formation. 

The parameters of the simulations are summarised in 
Table [1] Due to the computational load of the RT calcu- 
lation, we primarily use the N144L10 series for our main 
results. The other two series (N216L10 and N400L100) are 
used for the resolution tests, with a fewer number of simu- 
lated galaxies postprocessed with RT. The adopted cosmo- 
logical parameters are consistent with the WMAP results: 
Ho = 72kms"^Mpc~^ (h = 0.72), = 0.26, D.a = 0.74, 
= 0.044, (78 = 0.80, and Us = 0.96. 



2.2 Stellar Radiation Transfer 

To estimate /esc, we compute the stellar radiation transfer 
and the ionization structure of gas in each dark matter halo 
by post-processing the simulation output. First we set up 
a uniform grid around each dark matter halo with a grid 
cell size equal to the gravitational softening length of the 
simulation, and translate the SPH gas information into a 
gridded data. The grid typically have ~300^ cells for high- 
mass haloes (~ 10^^ Mq) and -70^ cells for low -mass haloes 

(~ wH'Iq). 

The RT scheme used in this paper is the Authentic 
Radiation Trans f er (A RT) method originally developed by 
iNakamoto et al.l (|200lf ). and the treatment is basically the 
same as in Y09. The performance of our RT code has already 
been repor t ed as part of the RT code comparison study by 
llliev et al] (|2006l ). and it can calculate the ion ization struc- 
ture precisely (see Figure 6 in llliev et aLlbOOd ). 

Usually the short-characteristic method is computation- 
ally cheaper than the long-characteristic method by a factor 
of Ng, where A^g is the grid size of RT calculation. In the 
short-characteristic method, the amount of calculation is re- 
duced by the interpolation of optical depth from the near- 
est grids. However the short-characteristic method suffers 
from an artificial photon diffusion effect. Our ART method 
is based on the long-characteristic method, and it is de- 
vised to reduce the calculation amount to a similar level 
as the short-characteristic method. Hence the ART method 
is suitable for the calculation of /esc for a large number of 
galaxies. As for the scatter ing of photons, w e employ the 
on-the-spot approximation (|Osterbrocklll989l ). in which the 
scattered photons are assumed to be absorbed immediately 
on the spot. 

In this work, the RT equation is solved along rays 
with uniform angular resolution from each source. The num- 
ber of ionizing photons emitted from the source stars is 
computed based on the theoretical spectral energy distribu- 
tion ( SEP) given by the p opulation synthesis code PEGASE 
v2.0 (jFioc fc Rocca|[l997^ . We take only the star particles 
that are younger than 10^ yrs as the sources of ionizing 
photons, and consider the effect of age and metallicty of 
the stellar population by interpolating the table generated 
from the result of PEGASE. We shoot the radiation rays 
in a radial fash ion from each star particles. We assume the 
ISalpeterl (|l955h initial mass function with the mass range of 
0.1 - 5OM0. 

Typically the postprocessing RT calculation takes 
about 100 hours for a large grid of 300^, and 1 hour for 
a small grid of 70^ on a single CPU. Our ART code is paral- 
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lelized by MPI, and has a high parallehzation efficiency. We 
process each star particle in parallel, and each CPU calcu- 
lates the radiation field from each star particle. In practice, 
we use ~ 1 — 128 CPUs simultaneously to process one halo 
with RT. 

2.3 Dust Attenuation 

We also include the effect of dust attenuation by distribut- 
ing the interstellar dust proportionally to th e metallicity, 
with a size distribution of nd{ad) oc a^^'^ jMathis et al.l 
119771 ). where ad is the radius of a dust grain. We adopt 
the dust grain size range of 0.1 — 1.0 pm as our fiducial 
model. The dust ma ss is calculated as ma = 0.01mg{Z/ZQ) 
(|Draine et al.l 120071 ') . where md, mg, and Z are the dust 
mass, gas mass, and metallicity in a grid. The density of 
a dust grain is assumed to be 3gcm~'^ like silicates. The 
dust opacity is given by dTdust ~ Q{u)na^ndd£, where 
Q{v), Od, rid, and d£ are the absorption efficiency factor, dust 
size, number density of dust grains, and path length, respec- 
tively. Since the assumed range of dust size is larger than the 
waveleng th of Lyman limit, w e assume Q{v) = 1 for ionizing 
photons (|Draine fc Ledl 19841 '). 

2.4 UV Background Radiation 

The baryonic gas in galaxies can be ionized by the UVB, 
and heated up to --^ 10'* K. It would be ideal to compute the 
RT of UVB as well as the stellar radiation, but in practice 
it is a very expensive calculation. 

Let us briefly explain why the UVB RT calculation is 
much heavier than the stellar radiation transfer. For stellar 
radiation, using on-the-spot approximation, we only calcu- 
late the RT along the angular rays between stars and grids 
with the dilution factor by the distance, whereas for UVB we 
have to calculate the RT of all angular rays. Therefore in the 
ART method, the number of rays that we have to calculate 
is A'^star xN4,xNeX Npi^th ~ A''star X Ng for stellar radiation, 
and for UVB. If > A''star, the calculation amount for 
UVB is larger than the stellar radiation. In practice, A^g is 
greater than A'atar by ~ 2 — 4 orders of magnitude. 

Due to this difficulty of UVB RT, usually a uniform 
UVB radiation field with an optically thin approximation 
is assumed across the simulation box as a simple approx- 
imation. Our fiducial simulation also in c ludes a uniform 
UV B with a modifie d iHaardt H Madaul (| 19961 ) spectrum 
(see iDave et al.l 1 19991 ). where the reionization takes place 
at z ~ 6 as sug gested by the quasar observations (e.g., 
iBecker et al ]|200ll) and stellar radiative transfer calculations 
(e.g.. lSokasian et al.|[20o3 ). 

However, the optically thin approximation is very crude, 
and the effects of different UVB has not been explored very 
much. Here we use following four N144L10 simulations with 
different treatment of UVB to examine the effects of UVB: 

(i) Fiducial: A uniform UVB radiation field with an 
optically thin approximation is assumed as stated above. 

(u) MHO. 5 (modified Haardt 0.5): The ISM is optically 
thin to the same UVB, however the intensity of UVB is 
reduced to the half of the Fiducial run. 



(ifi) OTUV (optically thick UV): The ISM is op- 
tically thin to the UVB in the lower density regions 
with Tin < O.Olpth = 6.34 x 10~^cm~^, but completely 
optically thick in higher density region (nn ^ O.Ol/Oth), 
where pth is the threshold density, above which star 
formation is allow e d. The value of pth was determined by 
IChoi fc Nagamin^ (|2009al ) based on the observed SF cut-off 
column density of the Kennicutt law. The OTUV method 
implicitly assumes that the UVB cannot penetrate into 
the high density regions by self-shielding. We find that 
this treatment reproduces the observed Hi column d ensity 
distribution function very well (|Nagamine et al.]|2010l ). and 
more detailed analyses using RT calculation supports this 
self-shielding density (Yajima et al. 2010, in preparation). 

(iv) no- UVB: UVB does not exist at all. 

We compute the ionization structure in each halo by 
solving the equation of ionization equilibrium as follows: 

(ruvB + Tstar) '^Hi + F*^ um ria = Qb nnii no, (2) 

where F^^g and T^^^^ are the photoionization rate by UVB 
and stellar radiation; F*^ is the collisional ionization rate; 
na,nm and nmi are the number density of free-electron, 
neutral and ionized hydrogen, respectively; qb is the total 
recombination coefficient to all bound excitation levels. The 
value of Fgtjj^ is estimated by the full RT calculation, but 
r"uvB computed with the optically thin approximation. 



3 RESULTS 

3.1 Ionization Structure 

Figure [T] shows the ionization structure of gas in a high- 
mass halo (Halo "A") and a low-mass halo (Halo "B") in 
the N144L10 Fiducial UVB run. 

The gas in Halo-A shows very complex and clumpy 
structure, going through continuous merging processes. The 
young star clusters are born in dense, neutral clumpy regions 
and irradiate the ambient ISM. However, the dense neutral 
gas clumps survive owing to high recombination rates. 

On the other hand, the Halo-B shows more or less spher- 
ical gas distribution before we process it with RT. Star clus- 
ters are born near the central high-density region. Once the 
halo is processed with RT, most of the low density gas on the 
right-hand-side is ionized by the UVB and stellar radiation. 
In particular, when the location of a young star cluster is 
slightly off-center, it can ionize one side of the halo preferen- 
tially, and creates a conical region of highly ionized region. 
The high density gas on the left-hand-side of the star cluster 
remains neutral, and the ionizing photons cannot escape to 
the left-hand-side. The value of angle-averaged /osc of each 
halo is basically determined by the covering fraction of these 
highly ionized region. 

3.2 Escape Fraction 

We estimate the average value of /osc for each dark matter 
halo as follows. For each light ray from each star particle in 
the halo, we count up the number of escaped ionizing pho- 
tons by integrating the transferred spectrum as a function 
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Figure 1. UpperiThe ionization structure of Halo A{Mh ~ 7 X IO^^Mq). lower:The ionization structure of Halo B{Mh ~ 1 X lO^^M©). 
Color indicates the neutral fraction of hydrogen in log scale. White points show the positions of young star clusters. 



of wavelength, and then divide it by the total number of 
intrinsically radiated ionizing photons. Then the values of 
/esc are averaged over all the rays coming out from the halo 
at the surface of the grid that was set up around the halo. 
Hereafter /esc denotes the angle-averaged value for each halo 
most of the time. 



3.2.1 Halo Mass and Redshift Dependence 

Figure[2]shows the /esc as a function of halo mass at z = 3 — 6 
for the N144L10 Fiducial UVB run. The solid triangles are 
the average values of /esc in each mass bin with l-cr error 
bars. We first take the average in the linear scale, and then 
take the logarithm for the data points (i.e., log(/esc)). At 
all redshifts, we find a clear qualitative trend that the mean 
/esc declines with increasing halo mass. 

Our results are similar to t he trend reported by 
iRazoumov fc Sommer-LarsenI (|2010l l. although our /esc val- 
ues are lower than theirs. Our results show the opposite 
trend to lGnedin et al.l (|2008l ). alt hough for high-mass galax- 
ies, our /esc is similar to theirs. In lGnedin et al.l (|2008l ). most 
of the simulated galaxies show a disk-like structure. In their 
scenario, for low mass galaxies, stars are born after the disk 



is formed, and most young stars are embedded deep inside 
the disk. As a result, most of the ionizing photons are ab- 
sorbed in the disk, and /esc is low. However in more massive 
galaxies, some star-clusters can form near the edge of dense 
disk, and they can be exposed by the mergers of galaxies. 
As a result of this efTect, they argued that more ionizing 
photons can escape from higher mass haloes. Therefore the 
geometry of simulated galaxies may cause the difference in 
the halo mass dependence of /esc. 

The galaxy sample size in iGnedin et al.l (|2008l ') and 
IRazoumov fc Sommer-LarsenI (|201Cll ) is ~ 10 — 20. and it is 
somewhat small to discuss the systematic tre nd of /esc as a 
function of halo mass. The spatial resolution of lGnedin et al.l 
(|2008l ). which is adaptively refined depending on the gas den- 
sity, is from ~ 17kpc to 260 pc. Although this resolution of 
maximum refinement level is better than that of ours, their 
RT schem e (OTVET) is co arser in estimating the ionization 
structure jlliev et al.ll2006l ). These differences in the accu- 
racy and resolution of fluid and RT calculations may have 
caused the difference in /esc. We will further discuss the pos- 
sible resolution effects in Section [5] 

Figure [3] shows t he redshift evolution of mean fgscin 
each halo mass bin. In|R azoumov fc Sommer- Larsenl(2010h . 
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Figure 2. Escape fraction as a function of lialo mass at z = 3 — 6 
for tlie N144L10 Fiducial UVB run. Different colors are used for 
different redshifts (red: 2: = 3, blue: z = b, green: z = 6). The 
triangles in the bottom right panel show the mean values in each 
mass bin with l-cr error bars. The data points with log /esc < —2.5 
arc shown at log /esc = — 2.5 for plotting purposes. 



the /esc of high-mass haloes with Mh > 10 Mq clearly de- 
creases with redshift (blue open circles) , and that of the low- 
mass haloes does not ch ange largely. On the other hand, our 
results and iGnedin et al. (20Q& ) indicate that /esc of high- 
mass haloes with Mu > 10^" does not change largely 
with redshift. For low-mass haloes with Mh < IO^^Mq, it 
seems that /esc is increasing slightly with decreasing redshift 
in our simulations. This might be due to the increasing cos- 
mic SFR density and increasing UVB intensity from z = 6 
to 2 = 3. Indeed, if we calculate the radiative transfer with- 
out the contribution of UVB in Eq. (2) for the Fiducial run 
at 2 = 3 with the same gas and stellar distribution, /esc 
decreases by ~ 10 — 20 per cent. In addition, the mass frac- 
tion of gas with lognn > 0.6 within haloes increases with 
increasing redshift, which leads to lower escape fraction due 
to higher recombination rate. 

Figure [4] shows the probability distribution function 
(PDF) of star particles as a function of /esc in haloes with 
Mh ^ 10" M0 (top panel) and Mh > lO^M© (bottom 
panel). The probability is defined by P(/esc) ~ A'^star(/esc ~ 
/esc + A/esc)/(A*'star,totaiA/osc), whcrc iVstar is the uumber 
of star particles that have the value of /esc, A^star, total is the 
total number of source star particles, and A/esc is the bin 
width. The figure shows that the lower mass haloes have a 
longer tail towards higher values of /esc. Since the ionization 
structure in low-mass haloes shows conical regions of highly 
ionized gas, ionizing photons can escape easily through these 
ionized cones, but not through other angular directions cov- 
ered by highly neutral gas. This allows for some star particles 
in lower mass haloes to have high /esc. On the other hand, 
the higher mass haloes show very complex and clumpy dis- 
tribution of highly neutral gas, therefore it is more difficult 
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Figure 3. Mean escape fraction in each mass bin of Fig. [2] as 
a function of redshift. Different symbols indicate different halo 
mass ranges (filled circles: lO* '^^ - IG^-^^M©, stars: lO^'^^ - 
W^-'^^Mq, filled triangles: lO^'^S _ iqIO-^SMq, open squares: 
1010-25 _ I0^0.75j^^^ crosses: IQiO-TS _ iQll-SSAfg, open trian- 
gles: 10"-25 - IO^^-'^^Mq, and open circle: 10^^ ''^ - 10^^-^^ Mq). 
The green open circles are for a relat ively massive halo (M tot^i = 
4 X W^^Mq at z = 3) examined bv iGnedin et al.l l|2008l '). which 
is to be compared with our open triangles. The red open circles 
ar e for a similarly massive halo (Mft, ~ 1O"M0) in the S33 run 
of lRazoumov fc Sommer^ Larsenllioioh . which is to be compared 
with our crosses. 



for the ionizing photons to escape, and there are no star par- 
ticles with /esc > 0.6. Thus the PDF for higher mass haloes 
is concentrated at /esc < 0.1. 

We also find that there is a large dispersion in /esc 
for the low-mass haloes with Mh ^ 10^^ Mq. This result 
may explain s o me of the recent observations . For example, 
IShaplev et"all (|2006l) and llwata et all (|2009l ) detected ion- 
izing radiation from high-2 galaxies, with a detection rate 
of about 10 per cent. The detected galaxies show extremely 
high /esc (~ 100 per cent). Our results do not show such 
high v alues of /esc, howe v er the /esc derived bv lShaplev et al.) 
(|2006l ) and llwata et all (|2009l ) are estimated from the flux 
ratio a t the Lyman limit and UV continuum. Recentlv llnoud 
(|2010t ) pointed out that the nebulae emission lines can boost 
the above flux ratio, leading to a very high /esc with an as- 
sumption that the /esc of nebular and stellar emission is a 
few tens of per cent. 

In our simulation sample, about 10 per cent show high 
/esc (> 0.4). These galaxies may corresponding to the re- 
cently observed o bjects with very high /esc- Furthermore, 
llwata et al.] (|2009l) showed that /esc decreases with increas- 
ing UV flux. In our simulations, SFR is positively correlated 
with halo mass, therefore our result of decreasi ng /esc with 
increasing halo mass is consistent with that of llwata et al.l 
(|2009l) . 
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Figure 4. Probability distribution function (PDF) of star parti- 
cles as a function of /esc at z = 3 in the N144L10 Fiducial UVB 
run. The top panel shows the PDF of all star particles in haloes 
with Mfi ^ IO^^Mq, and the bottom panel is for haloes with 
A'/fi > 10^^ Mq. Lower mass haloes contain more star particles 
with larger /esc- 



3.2.2 Dependence on the UVB Models 

Figure [5] shows /esc as a function of halo masses for different 
UVB model runs ai z — 3. Similarly to the Fiducial UVB 
model run, /esc decreases as the halo mass increases in all 
UVB models. Since the ionization fraction of gas should in- 
crease with the increasing UVB intensity, one would naively 
expect a higher /esc for the runs with stronger UVB in- 
tensities. However, when the UVB intensity increases, the 
stronger gas heating results in less efficient cooling of gas and 
hence less star formation. This reduces the number of ioniz- 
ing photons, and decreases the ionization fraction around the 
star-forming regions. Therefore these two competing effects 
counteract each other and self-regulate the ionization frac- 
tion, resulting in no significant differences in foac between 
different UVB runs, as shown in the bottom right panel of 
Figure \5\ 

We also considered the possibility that the sites of star 
formation might be different depending on the UVB inten- 
sity. Figure |6] shows the probability distribution function 
(PDF) of star particles as a function of relative distance 
"7?rei" between each star particle and the highest gas den- 
sity peak in the halo, normalized by the maximum radius of 

the halo (Rrel = Rsta.T / Rhalo) . To foCUS OU the scatter of /esc 

that we see for the lower mass haloes, here we selected only 
the low-mass haloes that include only one young star parti- 
cle. (There could be other star particles that are older in the 
same halo.) This figure shows that, when the UVB intensity 
becomes weaker, star-forming regions are more spread out 
to larger radii. Therefore the weaker UVB runs have more 
extended tails to larger radii compared to the Fiducial UVB 
model, progressively in the order of MHO. 5, OTUV, and no- 
UVB run, although the probability in the extended tail is 
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Figure 5. Escape fractions of ionizing photons (/esc) for the four 
UVB models at 2 = 3. Different colors indicate different UVB 
models (red: Fiducial, blue: MHO. 5, green: OTUV, and magenta: 
no-UVB). The triangles show the mean values in each mass bin 
with 1-(T error bars. The data points with log /esc < — 2.5 are set 
to log /esc = —2.5 for plotting purposes. There is not much dif- 
ferences between different UVB model runs due to self-regulation 
effect described in the text. 



very small. If a star cluster is farther away from the gas den- 
sity peak, the value of /esc would increase because the solid 
angle subtended by the high-density neutral clouds would 
become smaller. 



3. 2. 3 Origin of Scatter in /esc 

To further investigate the effect of star cluster locations on 
/esc , we plot /esc as a function of R^^i for each halo in differ- 
ent UVB models in Figure [T] Here again, we selected only 
the low-mass haloes that include only one young star parti- 
cle to focus on the scatter in /esc of the low-mass haloes. The 
gas near the density peak has high recombination rate, and 
therefore optically thick. The value of /esc strongly depends 
on the size of viewing angle to optically thick cloud from 
the source. As the source deviates from the central density 
peak, the viewing angle towards optically thick clouds de- 
creases, and on average, /esc increases with increasing Riei. 
Moreover the mean value of /esc does not depend on the 
UVB models very much, although the scatter is somewhat 
larger at larger i?rei values owing to the small sampling. At 
the lower values of Rrei, the scatter among different UVB 
models are smaller, because UVB cannot penetrate into the 
high-density cloud by the self-shielding effect. The results 
shown in Figures [S] and [7] suggest that the variation in i?rei 
is one of the key factors that determines the scatter in /esc 
for low-mass haloes. 

Furthermore fcsc may depend on the hydrogen number 
density at the location of star clusters, because the neutral 
hydrogen gas near the star clusters can effectively absorb 
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Figure 6. PDF of star particles as a function of ij^cl in different 
UVB models at 2: = 3, where R^d is the relative distance between 
each star particle and the highest gas density peak in the halo, 
normalized by the maximum radius of the halo. To focus on the 
scatter of /esc , here we used only the low-mass haloes that include 
only one young star particle. The result of the Fiducial UVB run 
is shown by the dashed histogram in other panels for comparison. 
The runs with a weaker UVB have more young star particles at 
larger distances from the density peak. 




Figure 7. Escape fraction of the low-mass haloes as a function of 
-Rrol. Similarly to Figure |6] only the low-mass haloes that include 
only one young star cluster are used. Different colors indicate 
different UVB models (red: Fiducial, blue: MH0.5, green: OTUV, 
and magenta: no- UVB). The triangles show the mean values in 
each i?rei bin with 1-cr error bars. The data points with log /esc < 
— 1.5 are set to log /esc = —1.5 for plotting purposes. 



ionizing photons. Figure |8] shows the /esc of low-mass haloes 
as a function of hydrogen number density tih at the location 
of star clusters in different UVB models. In this figure we 
use only the low-mass haloes that include one young star 
particle to focus on the scatter of /esc, similarly to Figures [6] 
and [7] We find that /esc steeply decreases with increasing 
local riH. When wh > 4cm~"^, more than 90 per cent of 
the emitted ionizing photons cannot escape from the galaxy. 
Once young stars born at such a high density region, most 
of ionizing photons are absorbed on the spot because of high 
recombination rate. Therefore we find that the variation of 
the local hydrogen number density also causes the large scat- 
ter in /esc for low-mass haloes. This result is consistent with 
that shown in Figure [71 because we expect lower riH for grid 
cells at greater -Rrei- 



3.2.4 Dependence on Other Physical Quantities 

We also examine the dependence of /esc on other physi- 
cal quantities, such as metallicity, SFR, and specific SFR 
(= SFR/stellar mass) in Figure |9l This figure includes all 
star-forming haloes, just like Figure [S] We find that /esc de- 
creases with increasing metallicity and SFR, and vice versa 
for specific SFR. These correlations are expected, because 
the metallicity and SFR are both positively correlated with 
galaxy stellar mass and halo mass. However, if the metallic- 
ity increases, the dust attenuation could have an extra effect 
on /esc , which we will discuss in the next section. 



4 DISCUSSION 

4.1 Dust Attenuation Effect 

Interstellar dust can decrease the /esc by absorbing the ion- 
izing photons. We evaluate the effect of dust attenuation on 
/esc by comparing the RT result with and without the dust 
treatment, as described in 12.31 Figure [TO] shows the values 
of /esc with and without the treatment of dust extinction, as 
a function of halo mass at z = 3. The reduction rate of /esc 
does not depend on the halo masses very much, and ranges 
from to 20 per cent, with an average of 14 per cent. 

If the dust-to-gas ratio is the same, the optical depth is 
roughly proportional to aj^, where ad is typical dust size. 
Note that dr = [Q-Ka^m^/ {4,-Ka^'^p/'i)\dl cx a~^^ , where 
and p are dust mass and dust density. If we change the dust 
grain size distribution to a smaller size range of 0.03—0.3 /im, 
the mean reduction rate increases to ~ 38 per cent at z = 3. 
The weak dependence on halo mass is perhaps because most 
of the star-forming regions are enriched close to the solar 
metallicity, irrespective of the host halo masses. 

Our reduction rate is somewhat smaller than that re- 
ported by Y09, which is probably owing to the difference 
in the volume occupied by the metal rich gas. In Y09, the 
model galaxy was an isolated system, and the ISM was glob- 
ally mixed by the shock from supernova explosion, because 
there were no further infall of pristine gas from intergalactic 
space. On the other hand, in the present work, pristine gas 
can accrete onto the haloes from intergalactic space, which 
reduces the volume fraction occupied by the metal rich gas 
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Figure 8. Escape fraction of the low-mass haloes as a function 
of hydrogen number density at the location of the star particle 
in different UVB models at 2 = 3. Similarly to Figure [B] only 
the low-mass haloes that include only one young star cluster arc 
used. Different colors indicate different UVB models (red: Fidu- 
cial, blue: MH0.5, green: OTUV, and magenta: no-UVB). The tri- 
angles show the mean values in each mass bin with 1-cr error bars. 
The data points with log /esc < —2.5 are set to log /esc = —2.5 
for plotting purposes. 



around star-forming regions. However, the reduction rate of 
/esc reported in Y09 are still in the 1-cr range of the present 

work. 

Sommer- LarsenI (|2010l ) and lGnedin et all 

/esc is 



[ Razoumov _ , _ _. . ^ 

l|2008l ) reported that the e ffect of dust a t tenua tion on /et_ 
very small. In particular, iGnedin et al.l (|2008h argued that 
in their simulations, only the ionizing photons from stars in 
the outer disk can escape from the halo, and the escaping 
photons pass through only the low density regions with low 
metallicity and low dust content. This may lead to smaller 
variations of /esc in their simulations compared to other sim- 
ulat ions that do not fully resolve the dis k structure. 

iRazoumov fc Sommer-LarsenI (|2010l ) achieved a spatial 
resolution of 0.1 — Ikpc in their SPH simulations using a 
zoom-resimulation technique, but it is not clear if they had a 
sufficient resolution to resolve the disk structure of galaxies. 
A lso, our SPH simulations have l ower resolution than those 
ofH azoumov fc Sommer-LarsenI (|20ld ). therefore we would 
expect a larger variation in /esc compared to theirs, and 
our results are consistent with this expectation. In our SPH 
simulations, the escaping photons may also pass through 
moderately high density regions with higher metal and dust 
content, with stronger effects of dust attenuation. 

Furthermore , our dust model is different from those 
of Gnedin et al.l (|2008l ) and IRazoumov fc Sommer-LarsenI 
(|2010D . They adopted dust extinction curves of Large and 
Small Magellanic Clouds, whereas we use the dust siz e dis- 
tribution of our Galaxy derived bv lMathis et al.l |l973) for a 
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Figure 9. Escape fraction as a function of metallicity (upper 
panel), SFR (lower left panel), and specific SFR (lower right 
panel) for all the star-forming galaxies in the N144L10 Fidu- 
cial UVB run at z = 3. The negative correlation between /esc 
and metallicity is mainly caused by the positive correlation be- 
tween halo mass and metallicity. It is not the metals that directly 
controls /esc- The data points with log /esc < —2.5 are set to 
log /esc = —2.5 for plotting purposes. 



silicate-type dust. Despite the fact that the effective absorp- 
tion cross section of our dust model is smaller than in their 
models, the scatter of /esc in our simulations may be larger 
than theirs owing to the differences in both the resolution 
and dust treatment. In the future, we plan to improve our 
dust model by including the treatments of formation and 
destruction of dust particles. 



4.2 Contribution to IGM Ionization 

As a result of our RT calculation presented in the earlier sec- 
tions, we are able to estimate the total number of ionizing 
photons that escape from all the star-forming galaxies in the 
entire simulation box. Our calculation includes haloes with 



Mh > W^-'^Mp; at 



3, and those with Mh > lO^Mp, at 



z = 6. Figure [TT] compares the comoving emission rate den- 
sity of ionizing photons Nion (i.e., the number of ionizing 
photons emitted per unit time and per unit volume) in our 
N144L10 Fiduci al UVB run to the required iVion to reion- 
ize the universe (jMadau et al.lll999l ). which is shown by the 
black solid curves for the clumping factors of C = 1,3, 10, 
and 30. The blue solid circles are for the intrinsically radi- 
ated photons from all star-forming galaxies in our simula- 
tion, and the red circles are for the escaped photons after 
the RT calculation. 

The clumping factor of IGM at z > 6 is still very un- 
certain, and the results from numerical simulations vary de- 
pending on the resolution and the treatment of physical pro- 
cess es such as star formation and radiation transfer. Ear- 
lier, iGnedin fc Ostrikeil (ll997^ suggested C = 30 at z ~ 6, 
however Iliev et al. I (|2007l ~ reported C = 10 us ng higher 
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Figure 10. Upper panel: Comparison of /esc with and without 
dust extinction in the Fiducial N144L10 UVB model at z = 3. 
Open black circles show /esc without dust extinction, and the 
filled red circles show fcsc with dust extinction. Magenta crosses 
show the case when the dust grain size range is changed to 
0.03 — 0.3/im. The triangles indicate the mean values in each mass 
bin with l-tr error bars. The data points with log /esc < —2.5 are 
set to log /esc = —2.5 for plotting purposes. Lower panel: Mean 
reduction rate of /esc for each mass bin when including dust ex- 
tinction. Different colors show different redshifts (red: 2 = 3, blue: 
z = 5 and green: 2 = 6). The dashed magenta result is when the 
dust grain size range is changed to 0.03 — 0.3/im. The error bars 
are l-cr. 



resolution simulatio ns with a RT treatment. More recently, 
iPawlik et all (|2009t ) reported C ~ 3 - 6 using a cosmological 
SPH simulation with an optically thin approximation. 

Our results show that Men of escaped ionizing photons 
after the RT calculation is greater than the required Mon 
to ionize the universe at z = 6 if C = 10, but below the 
required value if C = 30. Therefore our fiducial simulation 
suggests that the star-forming galaxies can ionize the IGM 
as long as C ^ 10. 

Our results on A^ion is higher than those derived by 
iBoIton fc Haehneltl (|2007l '). which was derived by using the 
results of cosmological SPH simulations and observational 
data of Lya forest. Although the error bars of their data 
points look very small, it represents only the dispersion of 
Lya opacity data. There are still significant uncertainties in 
the sp ectral shape and the mean free path of ionizing pho- 
tons in lBolton fc Haehneltl (|2007l ). In their calculation, they 
use the distance between Lyman limit systems as a mean free 
path of ionizing photons. However many low density H I gas 
clouds can decrease the mean free path, and hence increase 
A'ion. Together with the uncertainties in our simulation such 
as the resolution and the details of star formation and feed- 
back modelSj_the differences between our results and that of 



Figure 11. Comoving emission rate density of ionizing photons 
A'^ion as a function of redshift. Blue filled circles are for the intrin- 
sically radiated photons from all simulated galaxies, and the red 
filled circles are for the escaped ionizing photons after the RT cal- 
culation. Black solid lines show the required to reionize the 



universe, derived bv lMadau et al.l jiggg) for various clumping fac- 
tors of IGM. Black dotted line sh ows the contribution from QSOs 
estimated bv lMadau et al. The filled squares are Njq^ de- 

rived from the Ly-o opacity data of IGM bv lBolton fc Haehneltl 
1I2OO7I) . 



ing photons by the haloes with different masses, as shown 
in Figure [T2| For the intrinsically radiated photons (blue 
lines), there is no clear trend with the halo mass, and all the 
haloes contribute roughly equal number of ionizing photons. 
However, the figure shows that most of the escaped ionizing 
photons (red histograms) come from the lower mass haloes 
(^ 10^" M0). This is because in our simulations, higher mass 
haloes have lower /esc. 

Earlier in Sect ion | 3.2.1l we di s cussed the large variation 
of /esc derived by Shaplev et all ( 20061') . Based on a clus- 
tering ana lvsis. lAdelberger et al] ( 20051 ) reported that the 
sample in IShaplev et al.1 (|2006l ') are hosted by haloes with 
Mk > W^^Mq. When compared wit h our results i n Fig - 
ure 1121 it suggests that the sample of IShaplev et al.l (|2006l ) 
might not be tracing the bulk of ionizing sources at z = 3, 



iBolton fc Haehneltl (|2007l ') can be accounted for. 

We also study the fractional contribution to the ioniz- 



and only sampling the massive end of t he distributio n. 

In addition, the detected sample in llwata et al.l (|2009l ') 
is brig hter than ~ 2 7 mag in R band. If we use the equa- 
tion in iMadau et al.l (| 19981 ) and the relation between halo 
mass and SFR (M^ ~ SFR x 10"Mq), the limiting magni- 
tude of Iwata's sample corresponds to Mh ~ 1.5 x IO^^Mq 
(SFR ~ 1.5M0yr-i). Therefore they are also tracing only 
the massive end of the distribution, although one order o f 
magnitude deeper in halo mass than IShaplev et al.1 (|2006l ). 
Observations of fainter sources are needed to capture the 
entire ionizing radiation from lower mass haloes. 

In the SPH simulations used for this work, haloes with 
are not resolved well. For example, in the 
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Figure 12. Number fraction of ionizing pliotons contributed by 
the haloes of different masses. The blue lines are for the intrin- 
sically radiated ionizing photons, and the red histograms are for 
the escaped ionizing photons after the RT calculation. This fig- 
ure shows that most of the escaped ionizing photons mainly come 
from lower mass haloes. 



case of the Fiducial N144L10 run, the halo with Mh = 
2 X 10^ Mq consists of 100 dark matter particles. In the 
future, we plan to study fcsc of even lower mass haloes using 
higher resolution simulations. If these low-mass haloes are 
the primary sources of ionizing photons at z > 6, the H ll 
bubbles during the reionization epoch will be produced by 
numerous low mass haloes. In addition, the large dispersion 
in fesc that we found in this work suggests that the sizes of 
the H II bubbles may have a large variety in the early stages 
of reionization. 



4.3 Resolution Test 

In the present work, we find that the value of /esc depends 
strongly on the distribution of star particles with respect 
to the high density gas. However, the dumpiness of ISM 
around the star-forming regions could strongly depend on 
the resolution limit of simulations, therefore it is important 
to evaluate the resolution effects on /esc- In particular, the 
cosmological SPH simulations have a difficulty in resolving 
the clumpy structure of ISM when the number of SPH parti- 
cles is very small, and the values of fcsc for low-mass haloes 
may be more strongly affected by the limited resolution. 
When the dumpiness of ISM is very high, it may have both 
positive and negative effect on /esc: the positive effect is that 
the ionizing photons may be able to escape through the void 
regions more easily, however those photons soon may be ab- 
sorbed by the nearby high-density neutral clumps, which 
would be a negative effect. 

To study the resolution effect. Figure [T3l compares /esc 
in the three runs with different resolution, as described in 
Table[T] The lowest resolution run (N400L100) with a larger 
box size shows overall higher /esc than the other runs. The 
highest resolution run (N216L10) gives similar values of 
/esc to the fiducial resolution run (N144L10), but somewhat 
lower /esc for the lower mass haloes with Mh ^ lO^'^A/©. 
Given the this result, it is quite possible that our /esc results 
have not converged yet, and we might still be overestimating 
/esc for the low-mass haloes in the present work. Neverthe- 
less, it is true in all the runs that /esc has a large scatter 
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Figure 13. The resolution test on /esc using three runs with 
different resolution, as described in Tablo[T] The lower right panel 
compares the mean /esc in each mass bins for the three runs 
with 1-(T error bars. The points with log /esc < —2.5 are set to 
log /esc = —2.5 for plotting purposes. 



for the low-mass haloes, and that /esc decreases on average 
with increasing halo masses. 

In Figure 1131 for the N216L10 run, we had to confine 
the sample to the low-mass haloes with Mh < 10^^ Mq, 
because of the heavy computational load. In the N216L10 
run, we need to set up a large grid of > 500"^ for the high- 
mass halos with Mh ~ 10^^ Mq, and these grids take too 
long to process with RT when we want to process a large 
sample. We will tackle the systematic study of haloes with 
higher resolution simulations using the next generation of 
supercomputers . 

For the low-mass haloes, the /esc of N216L10 run is 
smaller by a factor of ~ 2 than in the N144L10 run. If we 
assume that /esc of all haloes becomes one half, the resulting 
Nion also becomes a half. Then the red circles in Figure [TT] 
would decrease by 0.3 dex for the N216L10 run, and the 
threshold dumpiness factor for IGM reionization changes to 
~ 10(3) at z = 3(6). 



5 SUMMARY 

We have performed three-dimensional radiation transfer cal- 
culations of stellar radiation for a large number of high- 
z star-forming galaxies in cosmological SPH simulations to 
explore the escape fraction of ionizing photons. Our major 
findings are as follows: 

• The value of /esc decreases steeply with increasing halo 
mass, irrespective of numerical resolution. 

• There is a large dispersion in /esc for low-mass haloes 
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with Mh ^ 10" M0. 

• The values of /esc do not vary much with redshift and 
different UVB models. 

• The average reduction rate of /esc owing to the dust 
attenuation effect is ~ 14% with a large dispersion. 

• The results of our Fiducial N144L10 run suggests that 
the star-forming galaxies can ionize the IGM at z = 3 — 6, 
if the clumping factor is C < 30 (10) at 2 = 3 (6). If we 
use the results of the N216L10 run, we roughly estimate 
that the above threshold values would change to C < 10 
(3) at 2; = 3 (6). Our results suggest that the star-forming 
galaxies become the main contributor of IGM ionization at 
3 < z < 6. 

• The low mass haloes with Mu < 10^" are the main 
ionizing sources of IGM in our simulations owing to their 
high /cBc- The fraction of escaped ionizing photons coming 
from the haloes with Mh < 10^° M© at « = 3 - 6 is 70 per 
cent for the Fiducial N144L10 run. 



As we summarised in Section [T] the current results 
on the escape fraction of ionizing photons are confus- 
ing, as different results are obtained from different sim- 
ulations. For example, iGnedin et al.l (|2008l ) argued that 
/esc increases with increasing halo mass in the range of 
Mh = 10^" - 10^^ A/©, and their values of /esc were 
mostly less than a few per cent, much smaller than the 
other published work. The trend found in our simulations 
(decreasing /e sc with increasing halo mass ) is si milar to 
that found bv lRazoumov &: Sommer-LarsenI (|2010t ). but we 
find lower /esc values despite of the fact that our simula- 
tions have lower resolution than their zoom-resimulations. 
Therefore the differences in /esc b etween our work and 
iRazoumov fc Sommer-LarsenI (|2010l ) cannot be explained 
simply by the resolutio n effect. 

We also note that IWise fc CenI (I2009D obtained much 
higher values of /esc (^0.4) than iGnedin et all l|2008l ') did, 
using the same AMR method, but for a different halo mass 
range. Considering these facts, the differences that we see 
now in the results of /esc may have to do more with the 
different treatment of radiation transfer and the UV back- 
ground radiation, rather than the resolution or numerical 
technique. However, more detailed comparisons are needed 
to make more definite statements. 

One of our main points is that the variation in /esc 
is caused by the different geometry of ISM distribution in 
the halo. Recentlv lAgertz et al.l (|2010l ) suggested that super- 
novae feedback and star formation efficie ncy can determin e 
the geometry of a disk galaxy (see also ISales et al.ll2010l ). 
Therefore the uncertainties in the treatment of star forma- 
tion, feedback, and radiation transfer are all important for 
the calculations of /esc , and we need to continue to improve 
these models through comparisons with future observations 
of high-2 galaxies. 

If we allow ourselves to speculate even further and com- 
bine all the current results mentioned above, it is possible 
that /esc has a peak at Mh ~ 10^ — 10^" Mq as a function of 
halo mass at 2: = 3 — 6. But this is highly speculative and 
by no means based on any definite physical arguments. 



The strength of our current work is the large sample 
size of galaxies that we processed with RT. Our simulations 
also adopt a new galactic wind model which produces more 
favorable results on the cosmic star form ation rate and the 
IGM statistics such as C iv mass density (|Ghoi fc Nagamin^ 
|2010D . Although it is difficult to address the exact effect 
of our wind model on /esc unless we process simulations 
with different wind models with RT calculation, our test 
calculations showed that the effect is not so strong, and our 
main conclusions of this paper should remain unchanged 
even if we modify the wind model slightly. We consider that 
the most significant results in the current work are the large 
scatter of /esc for the low mass haloes, and its decline with 
the increasing halo mass. The earlier works by other authors 
did not discuss the scatter among different haloes with a 
wide range of mass owing to their small sample size. 

At 2 ^ 7, even lower mass haloes with Mh < W^Mq 
may become the main sources of IGM ionization. However 
in such low-mass systems, the UV radiation of massive stars 
may influence the gas dynamics signiflcantly (|Wise fc CenI 
[2003). Simulations with higher resolution than presented in 
this paper are needed to follow the star formation in such 
low-mass systems, and we need to solve the hydrodynamics 
and radiation transfer simultaneously to examine the effect 
of radiative feedback. In the future, we plan to couple the 
RT with hydrodynamics and study the effects of radiative 
feedback by young stars and AGNs. 
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